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We investigate the hard-core Bose-Hubbard model with fully anisotropic long-range dipole-dipole 
interactions on square lattices. In our model, we assume that the dipole moments are oriented 
along a particular axis on the 2D plane. To treat this model exactly, we perform unbiased quantum 
Monte Carlo simulations using a hybrid algorithm of the worm algorithm and an 0{N) Monte 
Carlo method. We obtain the ground-state phase diagram that includes a superfluid phase and a 
striped solid phase at the particle density p = 1/2 in broad regions. The obtained phase diagram 
indicates that a supersolid state is unstable. We give the qualitative discussion of the reason from 
a perturbative treatment. Finite-temperature transitions to the phases are also investigated. For 
large dipole-dipole interactions, we observe a small p = 1/3 striped solid phase and incompressible 
regions adjacent to it. In spite of its incompressibility, the particle density increases as the chemical 
potential increases in the regions. This indicates the devil's staircase caused by the presence of 
numerous metastable states. 

PACS numbers: 03.75.Hh, 05.30.Jp, 67.85.-d 



Since the experimental realization of dipolar ^^Cr 
Bose-Einstein condensation (BEC)[3, the physics of the 
dipole-dipole interaction has received growing attention, 
because the anisotropic and long-range interactions in- 
troduce new phenomena. For instance, anisotropic col- 
lapse and complex dynamics such as d-wave symmetric 
explosion have been confirmed in dipolar ^^Cr BEC0- 
More recently, there are intensive experimental efforts to 
to produce systems of cold polar molecules, where the 
strength and direction of electric dipole moments are 
tunable by applying a static electric field From 
theoretical aspects, recent quantum Monte Carlo stud- 
ies revealed the ground-state phase diagram of the hard- 
core Bose-Hubbard model with dipole-dipole interactions 
on square lattices @ and triangular lattices 0- In these 
works, dipole moments are assumed to be oriented per- 
pendicular to the two-dimensional (2D) plane. Thus, 
the dipole-dipole interaction is the isotropic repulsive 
one; where r is the distance between two inter- 

acting particles. Remarkably, in this condition, the 
checkerboard supersolid on square lattices turns out to be 
stabilizedQ, although it is not only by nearest-neighbor 
and next-nearest-neighbor interactions!^ In contrast 
to the isotropic case, anisotropic interactions derived 
from alignment of dipole moments may produce other 
quantum phases. In the soft-core dipolar Bose-Hubbard 
model, the mean-field calculation based on a Gutzwiller 
ansatz predicted striped supersolid phases in 2D square 
lattices and the layered supersolid phase in 3D cubic 
lattices [lot. However, an unanswered question is whether 
anisotropic long-range dipole-dipole interactions also sta- 
bilize supersolid states in the hard-core bosonic case. 

Although the unbiased quantum Monte Carlo method 
is a powerful tool to investigate quantum many body sys- 
tems, there is a severe problem to perform simulations for 
long-range interacting systems such as dipolar systems. 
The difficulty is that the computational cost becomes 



0(-/V^), whereas it is 0{N) in short-range interacting sys- 
tems. Here, N is the system size. To overcome this diffi- 
culty which occurs in general Monte Carlo methods, Lui- 
jten and Blote proposed an Monte Carlo algorithm which 
enables simulations with 0(A^ log iV) costs even in the 
presence of long-range interactions |11|. Quite recently, 
Fukui and Todo developed more efficient algorithm to 
treat long-range interactions with 0{N) costs In our 
previous study [l^, we applied the 0(iV) method to the 
worm (directed loop) algorithm [14| - [la | which enables us 
to simulate bosonic lattice systems with remarkable ef- 
ficiency. Using this algorithm, we revealed the presence 
of two-types of peaks in the momentum distribution of 
the checkerboard supersolid state in the isotropic dipolar 
systems. In this paper, we investigate quantum phases 
of bosons with anisotropic dipole-dipole interactions on 
square lattices by exact quantum Monte Carlo calcula- 
tions. Our algorithm is a hybrid algorithm of the worm 
algorithm and the 0{N) method mentioned above. 

Systems of dipolar bosons in an optical lattice are de- 
scribed by the Bose-Hubbard model with the on-site in- 
teraction and dipole-dipole interactions [TtI ITsj. In the 
case of anisotropic dipole-dipole interactions on 2D sys- 
tem, the on-site interaction should be strong to prevent 
system collapse caused by the attractive part of dipole- 
dipole interactions. When we consider situations where 
the on-site repulsion is strong and the particle density is 
low, it is reasonable to treat the hard-core Bose-Hubbard 
model with dipole-dipole interactions for simplicity. The 
Hamiltonian that we consider is given by 



H 



h.c. 



{%]) 



where 



V 



|rjjp(dj • dj) - 3(dj • rij){dj ■ ri 



,(1) 



(2) 



2 



(a) 



a. 



-1.8 



-2 - 



■2.2 - 



-2.4 - 



-2.6 




(b) 



E 



y (attractive) 



(repulsive) 



0.62 



0.64 



0.66 



t/v 



FIG. 1: (Color online) (a) Ground-state phase diagram of 
liard-core bosons with anisotropic long-range dipole-dipole in- 
teractions on square lattices. The dipole moments are point- 
ing parallel to the y-axis. Error bars are drawn but most of 
them are smaller than the symbol size(here and the following 
figures). Dashed lines are schematic phase boundaries, (b) 
Schematic configuration of the stripe solid state at p = 1/2. 
Bosons are represented by simple circles. 
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FIG. 2: (Color online) Particle density p, compressibility 
K, structure factor S{n,0)/N , and superfluid stiffness ps 
as functions of the chemical potential p/V at {t/V,T/t) — 
(0.62,0.05). 



Here, 6|(6i) is the bosonic creation(annihilation) opera- 
tor on the site i, and is the particle number operator 
defined by rii = blbi. t, /i, and V{> 0) are the hopping 
parameter, the chemical potential, and the strength of 
the dipole-dipole interactions respectively. Vij = ri — rj 
is the relative coordination vector between the site i and 
j. di is the unit direction vector of dipole moment of 
the particle at the site i. To investigate the anisotropic 
nature of the above model, we focus on the case where 
electric (or magnetic) dipole moments are oriented in the 
y-axis on the 2D x — y plane by applying a static uni- 
form electric field E(pv magnetic field B). In this case, 
the direction of dipole moments is d = (0,1,0) regard- 
less of the sites. Thus, the dipole-dipole interaction Eq. 
^ becomes the form of Vij = (r^ — 3r^)F/r^, where r.y 
is the distance in the y-direction. Therefore, the system 
has attractive long-range interactions in the y-direction 
in contrast to repulsive ones in the x-direction. In our 
simulations, we treat the N = L x L square lattice sys- 
tems with the periodic boundary condition. The lattice 
spacing is set to unity. To eliminate the effect of cutoff 
in the long-range interactions, we employed the Ewald 
summation method [l^. 

Our main result is the ground-state phase dia- 
gram shown in Fig. [T] To obtain the ground- 
state properties, we calculate the particle density 
p = the compressibility k = dp/dp, = 

[((E»"0^) - (E^n^y]/iTN), the superfluid stiffness 
Ps = {W^)T/t, and the structure factor S{k) = 
1 J2i j e^^'''"'^ {{niTij)- (rii)'^) for a sufficiently low tem- 
perature T/t = 0.05. Here, (• • • ) indicates the thermal 
expectation value and W = {Wx,Wy) is the winding 
number vector in the world-line representation [23. In 
the ground-state phase diagram for t/V > 0.62, we found 



a superfluid (SF) phase and a Mott lobe at p = 1/2. The 
Mott lobe at p = 1/2 corresponds a striped solid (ST) 
phase which is characterized by finite value of ^(Tr, 0)/N, 
and vanishing n and ps- Fig. [2]shows plots of these quan- 
tities as functions of p/V at {t/V,T/t) = (0.62,0.05). 
Since all numerical observables show clear jumps, the 
boundaries between different phases are separated by 
first-order transitions. In 2D systems with isotropic dipo- 
lar interactions, there is a theoretical prediction that 
first-order transitions with a density change are forbidden 
due to the negative log-divergent surface tension between 
two phases |21|. In contrast, when the dipoles are point- 
ing in the 2D plane, the sign of the surface energy can be 
non-negative, and, therefore, first-order transitions are 
allowed j22). For smaller hopping parameters t/V < 0.61, 
we observed a small striped solid phase at p = 1/3 (not 
shown in Fig. IT]) and incompressible regions like devil's 
staircase (DS) [a [l^-HBj . In our simulations, we found no 
evidence of a striped supersolid phase. 

For small hopping parameters t, the absence of striped 
supersolids can be understood qualitatively by discussing 
the stability of the supersolid against domain wall 
formations [2g. Although we discuss the possibility of 
interstitial-based supersolid state below, the same argu- 
ment can be also applied to vacancy-based supersolid 
states because of the particle-hole symmetry in hard- 
core bosonic systems. In Fig. [3l we show a possible 
supersolid by delocalization of interstitials on the striped 
background ^27j ^] in Fig. (Ha), and a domain wall 
formed by doped particles in Fig. ^h). In both situ- 
ations, we assume that particles with density of ~ l/L 
are doped into the p = 1/2 striped solid state. We first 
consider the classical limit t = 0. When we focus on 
interactions between doped particles, we notice that the 
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FIG. 3: (Color online) (a) Possible supersolid by delocal- 
ization of interstitials (shaded circles) on the striped solid 
background (simple circles). Arrows indicates hopping pro- 
cesses, (b) A domain wall (dashed line) formed by doped 
particles(shaded circles). The wavy lines represent the attrac- 
tive nearest-neighbor interactions between doped particles. 
These interactions are the most strong interactions among 
the dipole-dipole interactions in the present case(V > 0). 



energetic cost of the domain wall formation is lower than 
that of delocalization. This is because the doped parti- 
cles gain large attractive energy by aligning in the attrac- 
tive direction, whereas, in the case of delocalization as in 
Fig. El^a), the interstitials are apart from each other, thus 
interacting weakly. Even in the presence of sufficiently 
small hopping parameter t, we can expect that the inter- 
stitials still form a domain wall, because of the energetic 
gain by the attractive interactions. Thus, with doping of 
infinitesimal particle density l/L, the supersolid state 
is unstable against the domain wall formation. When the 
hopping parameter t is increased, the situation becomes 
more complicated. This is because, when interstitials de- 
localize as in Fig. EJa), the kinetic energy gain is 0{t), 
while it is only 0{P) for the case of domain wall forma- 
tions as in Fig. Efb). This causes the possibility that the 
energy costs reverse for finite t. However, the absence of 
supersolid phase in our simulation results suggests that 
the energy cost of the supersolid state is still larger than 
domain wall formations even for finite t. 

We next investigate the finite-temperature transitions 
to the striped solid phase and the superfiuid phase respec- 
tively. To clarify the finite-temperature phase transition 
to the striped solid phase and its universality class, we 
calculate the structure factor S{k)/N and the Binder ra- 
tio g = l/2[3 — (77i^)/(to^)^], where m is the order param- 
eter defined by m = 1/Nj2i riiC^^''''' at the wave vector 
k — (7r,0). The results are shown in Fig. S) From the 
crossing point of the Binder ratio, we obtain the critical 
temperature as Tc/t = 0.0580(5) [see Fig. S^al)]. In the 
finite-size scaling analysis, we assumed the scaling forms 
[S{k)/N]L^>^^'' = firL^/") and g = HtL^/"), where / 
and h are scaling functions and t = {T — Tc)/Tc. Since 
the finite-temperature transition is related to a transla- 
tional symmetry breaking of Z2 in the repulsive direc- 
tion, the Ising-type universality class is expected. Using 
the above critical temperature and the critical exponents 
1/ = 1,[3 = 1/8 which belong to the 2D Ising universal- 
ity class, we successfully performed the finite-size scaling 
analysis as shown in Figs. Ill[a2) and|3l^b2). 

To discuss the finite-temperature transitions to the 
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FIG. 4: (Color online) Finite-temperature phase transition 
to the striped solid phase at p = 1/2. (al)Temperature de- 
pendence of the Binder ratio g and (a2) its finite-size scaling. 
(bl)Temperature dependence of the structure factor S{k)/N 
at fc = (tt, 0) and (b2) its finite-size scaling analysis. 
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FIG. 5: (Color online) Temperature dependence of the corre- 
lation ratio C{L/2, 0)/C(I//4, 0) for difi^erent system sizes. In 
the inset, finite-size scaling plots are shown. 



superfiuid phase, we calculate the correlation ratio 
C(L/2,0)/C(L/4,0), where C(r) = {hrbl). Fig. [5]shows 
the correlation ratio C(L/2, 0)/C(L/4, 0) as a function of 
the temperature. In this figure, we can confirm the merge 
of the data below a critical temperature, which is char- 
acteristic of the Kosterlitz-Thouless (KT) transitions [29l 
[30| . To estimate the critical temperature, we per- 
formed the finite-size scaling analysis for the KT tran- 
sitions. In this analysis, the scaling form is assumed to 
be C(L/2,0)/C(L/4,0) = /(i/exp[c/V(r - Tkt)/*]) 
where a constant value c and the critical ternperature 
Tkt should be determined simultaneously [sTi |32[ . The 
result of the scaling analysis is shown in the inset of 
Fig. [5l From this analysis, we estimated the unknown 
values as c = 1.17(27) and T^t/I = 0.334(11). We 
also confirmed a similar KT-like behavior for the ratio 
C(0, L/2)/C(0, L/4) and obtained the same critical tem- 
perature within error bars(not shown here). To clarify 
the anisotropy of the superfiuid phase, we show the mo- 
mentum distribution ri(fc) ~ '^l^^i jC{rij)e^^''''^' in 
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Fig. |6l It can be seen that anisotropics of the momen- 
tum distributions are strong as the hopping parameter 
becomes weak and, thus, the dipole-dipole interaction 
becomes relatively strong. A qualitatively similar behav- 
ior has been observed in ^^Cr BEC as anisotropic cloud 
obtained by the time-of- flight experiment [s^ . 
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FIG. 6: (Color online) (a)~(c) Momentum distributions of 
bosons n(fc) in superfluid states for different hopping param- 
eters at the linear system size L = 32. The chemical potential 
and the temperature are fixed at {^/V,T /t) — (—2.6,0.2). 
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FIG. 7: (Color online) (a) Particle density p and compress- 
ibility K as functions of the chemical potential fi/V for t/V — 
0.61. (b) Schematic configuration of the state at p = 1/3 
plateau. 

Finally, we discuss the results for t/V < 0.61. For 
t/V < 0.61, we found that the superfuhd phase disap- 
pears and incompressible regions appear instead of it. 
To show this, we plot the particle density p and the com- 
pressibility K as functions of the chemical potential fi/V 
in Fig. [Zlja). In the simulations, we employed the tem- 
perature annealing to prevent the simulations from being 
trapped in local minima. In contrast to the results for 
t/V > 0.62(see Fig. we observed vanishingly small 
(but finite) compressibility even between the empty re- 
gion and the p = 1/2 striped solid phase. In particular, a 
small plateau appears at p 1/3. As shown in the inset 



of Fig. [Tlja), the values of structure factor S{k)/N at 
fc = (2tt/3, 0) survives when the system size is increased. 
Therefore, the striped solid state at p ~ 1/3 has period- 
icity of 3 in the a;-axis. A schematic configuration of this 
state is shown in Fig. [7l[b). In the intermediate regions 
except for the plateau, in spite of the incompressibility, 
the particle density increases as the chemical potential 
increases. This characteristic feature of devil's staircases 
suggests the presence of the numerous metastable states. 
Unlike the case that the system has isotropic dipole- 
dipole interactions 0, it is naively expected that all or- 
derings in the devil's staircase can be explained by the 
stripe-type ones, because the dipole-dipole interactions 
along the y-axis are always attractive in the present case. 
In fact, from snapshots, we confirmed that, in the DS 
between p = 1/2 and 1/3, the configurations have mixed 
structm-es of the striped solid with periodicity 2 and 3 
like a fioating phase in the 2D ANNNI(anisotropic next- 
nearest-neighbor Ising) modeljlJl- Based on the broken 
symmetry in the DS and the SF, it is expected that the 
phase transition takes place. However, the precise anal- 
ysis of the phase boundary between the DS and the SF 
suffers from the strong system-size dependence. To dis- 
cuss the details, futher studies beyond numerical compu- 
tations are highly desiable. 

To summarize, we have investigated the hard-core Bose 
Hubbard model with the anisotropic dipole-dipole inter- 
actions by the unbiased quantum Monte Carlo calcula- 
tions. At low temperatures, there are several anisotropic 
phases, such as the superfluid phase, p = 1/2 striped solid 
phase, p = 1/3 striped sohd phase, and the devil's stair- 
case. In our simulation, supersolid phases have not been 
found. This is because doped particles into the solid pre- 
fer to form domain walls in the strong attractive direction 
instead of delocalization. Although we have treated hard- 
core bosonic systems with fully anisotropic dipole-dipole 
interactions, the striped supersolid might be stabilized 
in some other situations where the direction of dipolcs or 
the sign of V are changed and /or a finite on-site repulsion 
is included as suggested by previous studies [l0ll35|. 
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